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ABSTRACT 

The Fermi Large Area Telescope (LAT) reported the first definitive GeV detections of the binaries 
LS I +61° 303 and LS 5039 in the first year after its launch in June, 2008. These detections were unambiguous 
as a consequence of the reduced positional uncertainty and the detection of modulated 7-ray emission on the 
corresponding orbital periods. An analysis of new data from the LAT, comprising 30 months of observations, 
identifies a change in the 7-ray behavior of LS I +61°303. An increase in flux is detected in March 2009 and a 
steady decline in the orbital flux modulation is observed. Significant emission up to 30 GeV is detected by the 
LAT; prior datasets led to upper limits only. Contemporaneous Te V observations no longer detected the source, 
or found it -in one orbit- close to periastron, far from the phases at which the source previously appeared at 
TeV energies. The detailed numerical simulations and models that exist within the literature do not predict or 
explain many of these features now observed at GeV and TeV energies. New ideas and models are needed 
to fully explain and understand this behavior A detailed phase-resolved analysis of the spectral characteri- 
zation of LS I +61°303 in the GeV regime ascribes a power law with an exponential cutoff spectrum along 
each analyzed portion of the system's orbit. The on-source exposure of LS 5039 is also substantially increased 
with respect to our prior publication. In this case, whereas the general 7-ray properties remain consistent, the 
increased statistics of the current dataset allows for a deeper investigation of its orbital and spectral evolution. 
Subject headings: binaries: close - stars: variables: other - gamma rays: observations - X-rays: binaries - 
X-rays: individual (LS I +61°303) - X-rays: individual (LS 5039) 



L INTRODUCTION 

To date there are only a handful of X-ray binaries 
that have been detected at high (HE; 0.1-100 GeV) or 
very high-energies (VHE; >100 GeV): LS I +61°30 3 
(lAlbert et alJ l2006t lA cciari et al.' "2001 lAbdo et al.l l2009ah . 
LS5039 (lAharonian etal. 2005b; Abd o et al.l l2009b|) . PSR 
B 1259-63 (Aharonian et alj I2005at |Abdo et al. 2011»)^w 
X-3 (lAbdo et al., ,2009(1 Cvg X-1 (.Albert et al.. l2007t 
ISabatini et al.ll2O10l) . Recently, two new binaries were found: 
IFGL J1018. 6-5 856, with a period of 16.6 days found in 
the GeV regime (ICorbetet al.ll2011h and HESS J0632+05 7 
(iFalcone et al.ll201 Ih lOng et al. 11201 U iMariotti et a"ni201 Ih . 
for which a period of ^--^320 days was detected in X-rays (Bon- 
giorno et al. 2011). Of these sources only LS I +61°303, 
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LS 5039 and PSR B 1259-63 share the property of being bi- 
naries detected at both GeV and TeV energies. The other sys- 
tems have been unambiguously detected only in one band, 
either at GeV or at TeV, see, e.g., the case of Cyg X-3 in 
lAleksic et all (120101) . In the case of Cyg X-1, with the hint of 
TeV detection itself being at the level of 4 standard deviations 
(4(t), claims of detection at GeV energies by the Astrorive- 
latore Gamma a Immagini Leggero (AGILE) remain uncer- 
tain with concurrent Fermi Large Area Telescope (LAT) ob- 
servations ( iHill et al. 2010). It is yet uncertain whether these 
spectral energy distribution (SED) differences reflect an un- 
derlying distinct nature, or are just a variability signature in 
different bands. 

The nature of the binary compact object in LS I +61°303, 
LS 5039, HES S J0632+()57 and IFGL J1018.6-5856 is as yet 
undetermined dHill et alj |2010). Both neutron star (e.g. PSR 
B 1259-63) and probable black hole (e.g. Cyg X-3) binary 
systems have been detected at GeV energies and so both types 
of compact object are viable in the undetermined systems. 
Recently the Burst Alert Telescope (BAT) onboard Swift re- 
ported a magnetar-like event which may have emanated from 
LS I +61°303 ( jBarth elm y et al. 2008; Torr es et al. 20jJ). If 
true this would be the first magnetar found in a binary system. 

The early LAT reports of GeV emission from LS 5039 and 
LS I +61° 303 were based upo n 6-9 months of survey ob- 
servations (lAbdo et al.ll2009allbh . Both sources were detected 
at high significance and were unambiguously identified with 
the binaries by their flux modulation at the co rresponding or- 
bital periods, 26.4960 days for LS I +61°303 (lGregorvl2002l) 
and 3.90603 days for LS5039 (ICasares et al.l 120051) . The 
modulation patterns were roughly consistent with expecta- 
tions from inverse Compton scattering plus 7-7 absorption 
models, and were anti-correlated in phase with pre-existing 
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TeV m easurements (e.g.. [Albert et al.ll2009l : lAharonian et al.l 
120061) . The anti-correlation of GeV-TeV fluxes is in fact a 
generic feature of these models, where the GeV emission is 
enhanced (reduced) when the highly relativistic electrons seen 
by the obser ver enco unter the seed photons head-o n (rear- 
on); (e.g.. seelBedn arek 2007; Sierpowska-Bartosi k & Torres! 
120071; iDubuset al. 2008; Khangulyan et al. 200^). Fermi- 
LAT measurements provided a confirmation of these predic- 
tions. 

The spectra of both sources were best modeled with ex- 
ponential cutoffs in their high-energy spectra, at least along 
part of the orbit. Specifically, an exponential cutoff was 
statistically a better fit to the SED compared with a pure 
power law at phases surrounding the superior conjunction 
(SUPC) of LS 5039and in the orbitally averaged spectrum of 
LS I H-61°303. Statistical limitations of the data prevented the 
determination or the ruling out of an exponential cutoff in any 
part of the orbit of LS I H-61°303 or in the inferior conjunction 
(INFC) of LS 5039. The spectral energy distributions with the 
exponential cutoffs that were reported were remini scent of the 
many p ulsars the LAT has detected since launch dAbdo et alj 
l2009bh . although this was far from a proof of their pulsar na- 
ture. To date no pulsations have been found at GeV energies, 
or at any other wavelengths, despite deep dedicated searches 
(see e.g., Rea et al. 2010, 201 1). 

Since Fermi was launched, both the Major Atmospheric 
Gamma-ray Imaging Cherenkov Telescopes (MAGIC) and 
the Very Energetic Radiation Imaging Telescope Array (VER- 
ITAS) have performed observations of LS I H-61°303. No TeV 
detection was reported after October 2008, un til the source 
unexpectedly reappear ed, once, at periastron jAcciari et alj 
1201 UlOng et al. I2010t). At the sa me time, a hard X-ray multi- 
year analvsis dzhanp et al. Il2010l) and a long-term X-ray cam- 
paign on LS I H-61°303 using the Rossi X-ray Timing Explorer 
(RXTE) has been cond ucted covering the whole extent o f the 
LAT observations (see [Torres eTaniioiOt ILi et al.ll201 Ih . In 
addition, simultaneous and archival data from long-term mon- 
itoring of radio and Ha emission is available for comparison 
in a multi-wavelength context. In this work we present the 
results of the analysis of 30 months of LAT survey observa- 
tions of both LS I H-61°303 and LS 5039. We investigate the 
long-term flux variations of the sources, as well as variations 
in the amplitude of their orbital flux modulation, and we ex- 
plore the possible spectral variability for both systems, finally 
putting and interpreting these observations in the context of 
the source behavior at other frequencies. 

2. OBSERVATIONS AND DATA REDUCTION 

The Fermi-LALT is an electron-positron pair production tele- 
scope, featuring solid state silicon trackers and cesium io- 
dide calori meters, sensitive to photons from ~'20MeV to 
>300GeV (lAtwood et alJl2009l) . It has a large -2.4 sr field 
of view (at 1 GeV) and an effective area of ~8000 cm^ for 
>lGeV. 

2.1. Dataset 

The Fermi survey mode operations began on 2008 August 
4; in this mode, the observatory is rocked north and south 
on alternate orbits to provide a more uniform coverage, so 
that every part of the sky is usually observed for '^30 minutes 
every —3 hours. Therefore, the two sources of interest were 
monitored regularly without significant breaks, allowing us to 
draw a complete picture of their behavior in 7-rays over the 



last two years. The dataset used for this analysis spans 2008 
August 4 through 20 11 January 24. 

The data were reduced and analyzed using the Fermi Sci- 
ence Tools v9r20 packaged The standard onboard fil- 
tering, event reconstruction, and classification were applied 
to the data (iAtwood e t al. 20ilS)- The high-quality "diffuse" 
event class was used together with the "Pass 6 v3 Diffuse" in- 
strument response functions (IRFs). Time periods when the 
target source was observed at a zenith angle greater than 105° 
were excluded to limit contamination from Earth limb pho- 
tons. Where required in the analysis, models for the Galactic 
diffuse emission (gll_iem_v02.fit) and isotropic backgrounds 
{isotropic_iem_v02.txt) were usecPl. 

2.2. Spectral analysis methods 

The binned maximum-likelihood method of gtlike, in- 
cluded in the ScienceTools, was used to determine the in- 
tensities and spectral parameters presented in this paper We 
used all photons with energy >100MeV in a circular region 
of interest (ROI) of 10° radius centered at the position of 
LS I H-61°303 and LS 50 39, respectively. F or source mod- 
eling, the IFGL catalog (lAbdo et al.ll2010ah . derived from 
11 months of survey d ata, and the first Fermi pulsar cata- 
log (lAbdo et alJl20T0bl) were used; all sources within 15° of 
the ROI center were included. The energy spectra of point 
sources included in the catalog within our ROI are modeled 
by a simple power law. 
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with the exception of known 7-ray pulsars, which were mod- 
eled by power laws with exponential cutoffs described by; 
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The spectral parameters were fixed to the catalog values ex- 
cept for the sources within 3 degrees of the candidate loca- 
tion. For these latter sources, the flux normalization was left 
free. All of the spectral parameters of the two subject bina- 
ries were left free for the fit. Source detection significance is 
determined using the Test Statistic value, TS = -2 ln(Lo/i'i) 
which compares the likelihood ratio of models including, e.g., 
an ad ditional source, wit h the null-hypothesis of background 
only jMattox et alj|1996l) . 

To estimate the systematic errors, which are mainly caused 
by uncertainties in the effective area and energy response of 
the LAT as well as background contamination, we use the 
so-called "bracketing" IRFs. These are IRFs with effective 
areas that bracket those of our nominal IRF above and be- 
low by linearly connecting differences of (10%, 5%, 20%) at 
log(E/MeV) of (2, 2.75, 4) respectively. 

2.3. Timing analysis methods 

Lightcurves are extracted using aperture photometry, taking 
an aperture radius of 1° and using the gtbin tool. The ex- 
posure correction is performed with the tool gtexposure 
assuming the spectral shape of the source to be a power law 

See the Fermi Space Science Center (FSSC) website for details of the 
Science Tools: http://fermi.gsfc.nasa.gov/ssc/data/analysis/ 

*^ D escriptions of the models are available from the FSSC: 
|http://fermi. gsfc.nasa.gov/ssc | 
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Fig. 1 . — Left: Spectral energy distribution of LS 5039 as observed by the Fermi-hAI (red). The gray open and filled circles are H.E.S.S. spectra in INFC and 
SUPC, respectively ( Aha ronian et alj2Q06i) and are not simultaneous with the GeV measurements. Right: Fermi-LAT and H.E.S.S. spectra of LS 5039 in SUPC 
(blue) and INFC (red). 



with an exponential cutoff (see §§ l3.1ll4n i. These lightcurves 
are not background subtracted. The folded lightcurves shown 
in the subsequent sections are derived by performing gt 1 i k e 
fits for each phase bin. Therefore, all of them are effectively 
background subtracted. We check that both methods for gen- 
erating lightcurves, aperture photometry and gt 1 i ke fits, are 
consistent with each other when the former lightcurves are 
background-subtracted too. 

The primary method of timing analysis employed searches 
for periodic modulation by calculating the weight ed pe- 
riodogram of the ligh tcurve CLomb, .1976;; Scargld IT9821: 
ICorbet & DuboisI 120071) . The lightcurve is constructed by 
summing, for each photon, the estimated probability that the 
photon came from the source of interest. The probability 
will be both spatially and spectrally dependent. Because this 
technique allows for the correct weighting of each photon it 
intrinsically improves the signal-to-noise and allows the use 
of a larger aperture. This method has successfully been ap- 
plied to increasing the LAT sensitivity for the detection of 
pulsars (lKenil20l"lT) . However, in the basic form of this tech- 
nique, the weight for any particular energy/position is fixed. 
This means that changes in source brightness will not be re- 
flected in the weights and can result in incorrect probabilities. 
The calculation of probabilities was performed using the tool 
gt s r cpr ob and the same source model file derived from the 
IFGL catalog and used in the spectral analysis. Since the ex- 
posure of the time bins was variable, the contribution of each 
time bin to the power spectrum was weighted based on its rel- 
at ive exposure. Period error s are calculated using the method 
of lHorne & BaUunasI (119861) . 

3. LS 5039 RESULTS 

LS 5039 is located in a complicated region toward the in- 
ner Galaxy with high Galactic diffuse emission and many sur- 
rounding 7-ray sources. In particular, the LAT detected a 
bright (8.7 x 10"^ ph cm"^ s"' above 100 MeV) 7-ray pul- 



sar, PSR J1826-1256, - 2° away from LS 5039. Following 
the ana lysis performed in the earlier LAT paper (lAbdo et alJ 
l2009bh we discarded events whose arrival times correspond 
to the peaks of the pulsar cycle of PSR J 1826- 1256 in or- 
der to minimize the contamination from the pulsar. The ex- 
cluded pulse phase of PSR J1826-1256 i s 0.05 < < !) ,, < .2 
and 0.625 < 0p < 0.775 (see Fig. 37 in iRay et all mTH) ). 
which results in a loss of 30% exposure on LS 5039. To ac- 
count for the loss, a scaling factor of 1/0.7 is multiplied to 
fluxes obtained with maximum likelihood fits. 

3.1. Orbitally averaged spectrum 

The orbitally-averaged spectrum of LS 5039 was initially 
investigated by fitting a power law and a power law with an 
exponential cutoff to the data. We compare two m odels uti- 
lizing the likelihood ratio test dMattox et al.lll996l) . i.e. for 
the ratio 2 x Alog(Likelihood) we assume a x^-distribution to 
calculate the probabilities taking into account the correspond- 
ing degrees of freedom (Eadie et al. 1998). In this case, the 
significance of a spectral cutoff was assessed by comparing 
the likelihood ratio between the power law and cutoff power 
law cases, which is -2 ln(LpL/Lcutoff) = 94.9, where LpL and 
i'cutoff are the likelihood value obtained for the spectral fits 
with a power law and a cutoff power law, respectively. This 
indicates that the simple power law model is rejected at the 
9.1(7 level in favor of the cutoff power law. The best-fit param- 
eters for the cutoff power law model are T = 2.06 ± 0.06stat ± 
0. 1 Isyst and iicutoff = 2.2 ± 0.3stat ± 0.5syst GeV with a flux of 
F>ioo MeV = (6.1 ±0.3stat±2.1syst) X 10"'' cm"" s"' integrated 
above 100 MeV. Using a cutoff power-law spectral model, the 
maximum likelihood fit yields a test statistic of TS = 1623 for 
the LS 5039 detection; equivalent to ~ 40 cr. We also tried 
a broken power-law spectral model for LS 5039 in addition 
to an exponentially cutoff power law. We found the broken 
power law gives lower TS values than the exponentially cut- 
off power-law case. 
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Spectral points in each energy band were obtained by di- 
viding the dataset into separate energy bins and performing 
maximum likelihood fits for each of them. The resulting spec- 
tral energy distribution (SED) is plotted in Figure [T] together 
with the best-fit cutoff power law model. Interestingly, the 
SED shows significantly higher flux (one spectral data point) 
at > 10 GeV than the expected flux from the best-fit cutoff 
power law, possibly suggesting another component at high en- 
ergies. 

One idea explored for these sources (especially, for 
LS5039, see e.g., Torres 2010) is that the 7-ray emission 
could be understood as having two components: one would 
be the magnetospheric GeV emission from a putative pulsar 
and the other from the inter-wind region or from the pulsar 
wind zone. The latter would be unpulsed and would vary with 
the orbital phase, the former would be steady and pulsed. The 
current data for LS 5039 would indeed allow for this possibil- 
ity, especially because of the possible high energy component 
found in the GeV spectrum. 

To test the significance of the additional component, we 
added to the model a power-law source at the location of 
LS 5039 in addition to the cutoff power law source. Model 
A is a power law with a cutoff and model B a power law with 
a cutoff plus an additional power law. Thus, model B has two 
more free parameters compared to model A. According to the 
likelihood ratio test, the probability of incorrectly rejecting 
model A is 6.1 x 10"^ (4.5(t). 

The best-fit parameters for the putative additional com- 
ponent are F = 1.6 ± 0.4stat ± 0.3syst and f>io oev = (1-6 ± 
0.6smt ± 0.7,y,,) X 10"'" cm-- s"'. The addition of the high- 
energy component slightly affects the parameters of the cutoff 
power law. The best-fit parameters for the latter are F = 2.02 ± 
0.06stat±0.10syst and£'cutoff = 2.0±0.3stat±0.4syst GeV with 

a flux of /^>iooMeV = (6.0 ± 0.3stat ± l-9syst X 10"^ cm"^ s"' 
integrated above 100 MeV 
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Fig. 2. — Light curve of LS 5039 in time bins of six orbital cycles between 
lOOMeV-300 GeV. One orbit is 3.90532 days long. 



3.2. Phase-resolved analysis 

Following the H.E.S.S. analysis bv lAharonian et af] ( 120061) 
as well as the previous one, the whole dataset was divided into 
two orbital intervals: superior conjunction (SUPC; cj) < 0.45 
and 0.9 < 0) and inferior conjunction (INFC; 0.45 < < 0.9). 
The SUPC and INFC data were analyzed in the same way as 
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Fig. 3. — Lomb-Scargle power spectrum of the whole GeV data set on 
LS5039. 



the orbitally averaged data. Being consistent with our pre- 
vious paper, the power-law assumption for the SUPC spec- 
trum is rejected with -2 ln(LpL /Lcutoff) = 81.2, or at a rejec- 
tion significance of ^ 9cr. The best-fit parameters are F = 
2.07 ± 0.07stat ± O.OSsyst, ^cutoff = 1 .9 ± 0.3stat ± 0.3syst GeV 
and F>ioo MeV = (7.8 ±0.4s,at± 1.9syst) x 10"^ cm"^ i"'. 

Although a single power law was not rejected for INFC in 
our previous analysis using 10 months of data (lAbdo et alJ 
l2009b), a cutoff power law is preferred also for INFC with 
the present dataset. The likelihood ratio for the INFC data 
is -2 ln(LpL/Lcutoff) = 21.7, which corresponds to 4.7(t. The 
parameters for the INFC spectrum are F = 1.99 ± 0.1 3stat ± 
0.07syst, ^cutoff = 2.6 ± 0.7stat ± 0.9syst GeV, and F>ioo MeV = 
(3.9±0.4,tat± l.Ssyst) X 10"^ cm-2 s"'. Therefore, the SUPC 
and INFC spectral shapes are completely consistent with one 
another within the errors. The only difference is the normal- 
ization and hence the total flux. On the other hand, the spec- 
trum for INFC (red points in the right panel of Figure[T]) seems 
to exhibit additional structure below 1 GeV. The limited statis- 
tics and the large contribution of diffuse emission at low ener- 
gies, however, prevents solid conclusions on whether a more 
complicated fit (e.g. a double broken power law or a broken 
power law with a cutoff) would be preferred. 

We also searched for emission from the high-energy com- 
ponent in the SUPC and INFC spectra. However, the TS of the 
additional components compared with a power law with ex- 
ponential cutoff are only 13.6 and 10.9 for SUPC and INFC, 
respectively, and do not confirm a second spectral component. 
The SUPC and INFC SEDs were obtained using the same 
method as the orbitally averaged spectra and are plotted in 
the right panel of Figure[T] 

3.3. Light curve 

Figure |2] shows the lightcurve for LS 5039 over 30 months 
derived by performing gtlike fits on time bins which con- 
tain 6 orbital cycles each. The lightcurve for LS 5039 does 
not show any significant flux changes. Constructing the pe- 
riodogram of the weighted photon lightcurve yields a signifi- 
cant detection of a periodicity at 3.90532 ± 0.0008 days. This 
is consistent with the known orbital period of LS 5039. The 
Lomb-Scargle power spectrum of LS 5039 is shown in Figure 
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[3] The stability of the orbital modulation was investigated and 
no significant variation in the modulation fraction as a func- 
tion of time was found. 

4. LS I +61°303 RESULTS 
4. 1 . Orbitally-averaged spectral analysis 

We have derived the spectrum of orbitally-averaged LAT 
data, i.e., without any selection criteria (cuts) concerning the 
orbital phase, for the LS I H-61°303 system. The spectral 
points and corresponding best-fit using the updated dataset de- 
scribed in §|2]are shown in Figure |4] together with previously 
derived results from the LAT and TeV observations. Two 
sets of TeV data are plotted: we show the non-simultaneous 
data points obtained by the Cherenkov telescope experiments 
MAGIC and VERITAS. (These data correspond to phases 
around 0.6-0.7 and represent several orbits observed in the 
period 2006-2008, before Fermi was launched). Addition- 
ally, we show the latest measurements performed by VERI- 
TAS, which established a 99% C.L. upper limitH The new 
VERITAS upper limit spans several orbits during which, si- 
multaneously with our LAT data, no detection was achieved. 
The LAT data along the whole orbit are still best described by 
a power law with an exponential cutoff. The TS value for a 
source emitting 7-rays at the position of LS I H-61°303 with an 
SED described by a power law with an exponential cutoff is 
highly significant. The relative TS value comparing a fit with 
a power law and a fit with a power law plus an exponential 
cutoff clearly favors the latter, at the ~20(t level. The pho- 
ton index found is F = 2.07±0.02stat±0.09syst; the flux above 
100 MeV is (0.95±0.01stat±0.07syst) x lO'^phcm-^s"', and 
the cutoff energy is 3.9 ±0.2stat±0.7systGeV. Results for the 
obtained TS values for each fit to different datasets are listed 
in Table[T]and all fit parameters obtained for the exponentially 
cutoff power law models are listed in Table |2] 

Figured shows that the data point at 30GeV deviates from 
the model by more than 3a (power law with cutoff, red line). 
Although in our representation it is only one point, it is in it- 
self significant, with a TS value of 67 corresponding to ^%a. 
Therefore, and similarly to the case of LS 5039, with the 
caveat of having only one point determined in the SED be- 
yond the results of the fitted spectral model, we investigate 
the possible presence of a second component at high ener- 
gies. As in the case for LS 5039, we use the likelihood ratio 
test to compare two models: Model A is a power law with 
a cutoff and model B a power law with a cutoff plus an ad- 
ditional power law. According to this test, the probability 
of incorrectly rejecting model A is 5.7 x 10"'^ (7.8(t). The 
TS value for this extra power law component as a whole is 
172, larger than in the case of LS 5039, and its parameters are 
F = 2.5 ±0.3stat and F>io cev = (3.5 ±0.6stat) x 10"'° cm'^ s"'. 
The addition of the high-energy component affects the pa- 
rameters of the cutoff power law. The best-fit parameters for 
the latter, when including the former, are F = 2.00 ± 0.03stat 
and ^cutoff = 2.7 ± 0.3stat GeV with a flux of F>ioo MeV = 
(0.88 iO.OSstat) X 10"^ cm-2 s"'. Compared with the cutoff 
energy of 3.9 ±0.2 GeV we obtain after fitting only a power 
law with an exponential cutoff to the data, the cutoff energy 
decreases when the additional high-energy power-law compo- 
nent is introduced. 

We derive this differential upper limit by using the VERIT AS-reported 
integral flux upper limit for phases 0.6-0.7 JAcciari et alj|201 H) assuming a 
differential spectral slope of 2.6. 



10-^2 



\ — M M M 1 — III 

A 


• Fermi, 30 months 
A Fermi, 8 months 

• MAGIC, 2005/06 

VERITAS, 2006/07 
VERITAS, 2008/10 


* V 

\ 


J 
















_ T 

iiiiT* 






\ 








1 


III 


^ 1 1 


1 







10^ 10" 10^ 10^ 10^ 

Energy [MeV] 



Fig. 4. — The overall 30-months LS I +61°303 spectrum (red) in compari- 
son with the earlier published one (black) over 8 months is shown. TeV data 
points taken by MAGIC and VERITAS are shown in gray. They are not si- 
multaneously taken with the GeV data. The VERITAS upper limit (in green) 
is. 
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Fig. 5. — Lightcurve of LS I +61°303. Each point represents one orbit of 
26.496 days, whereas the black solid line represents the 3-bin smoothed light 
curve. The black dashed line marks the moment of the flux change in March 
2009. 

An alternative model to accommodate the deviating high- 
energy point in the spectrum is a fit with a broken power-law. 
This is further discussed in Section IT?! 

4.2. Lightcurve 

In Figure|5]the lightcurve for LS I -i-61°303 over 30 months 
is shown using orbital time bins. The black dashed line rep- 
resents the point in time when a flux change occurred for 
LS I H-61°303; this ocurred in March 2009. In Figure |6] we 
present the power spectrum of LS I -1-6 1 °303 derived from the 
total weighted photon lightcurve. The power spectrum clearly 
detects the orbital flux modulation with a period of 26.71 ± 
0.05 days. This is consistent with the known orbital period of 
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TABLE 1 

TS VALUES FOR LS I +61° 303 FOR DIFFERENT SPECTRAL SHAPES (SEE SECTION|0]fOR DETAILS) 



Data set 


Power law + exponential cutoff 


Power law 


Broken power law 




TS 


TS 


TS 


30 months of data 


23995 


23475 


23970 


Data before March 2009 


3404 


3314 


3415 


Data after March 2009 


20714 


20283 


20699 


Inferior conjunction (geometrically) 


12548 


12326 


12512 


Superior conjunction (geometrically) 


I171I 


11422 


11700 


Inferior conjunction (angle cut) 


6670 


6562 


6665 


Superior conjunction (angle cut) 


6083 


5986 


6063 


Periastron 


11656 


11450 


11636 


Apastron 


12377 


12059 


1 236 1 



TABLE 2 

PARAMETERS FOR LS I +61°303 FROM SPECTRAL FITTING WITH POWER LAW WITH EXPONENTIAL CUTOFF (SEE SECTIOn [431 fOR DETAILS) 



Data set 


Photon index V 


Cutoff energy 


Flux >IOOMeV 








[GeV] 


[xlO^^h cm"- s"'] 


First 8 months of data 


2.21 ±0.04|s,„ 


± 0.06,,,, 


6.3 ± 1. 1,,,., 


±0.4,y„ 


0.82 ±0.03,,., 


±0.07,y„ 


30 months of data 


2.07 ± 0.02,iai 


±0.09.ya 


3.9 ±0.2,,., 


±0.7,y„ 


0.95 ±0.01,,., 


± 0.07,y„ 


Data before March 2009 


2.08±0.04„ai 


± 0.09sys, 


4.0 ±0.6,,., 


±0.7,y„ 


0.75 ± 0.03,,., 


±0.07,y„ 


Data after March 2009 


2.07±0.02.„a, 


±0.09.ys. 


3.9 ±0.3,,., 


±0.7,y„ 


1.00 ±0.01,,., 


±0.07,y„ 


Inferior conjunction (geometrically) 


2.14±0.02aa, 


±0.09sys. 


4.0 ±0.4,,., 


±0.7,y„ 


1.07 ±0.02,,., 


±0.07,y„ 


Superior conjunction (geometrically) 


2.02 ±0.03,,., 


±0.09.ys. 


3.9 ±0.3,,., 


±0.7,y„ 


0.85 ±0.02,,., 


±0.07,y„ 


Inferior conjunction (angle cut) 


2.17±0.03.„a, 


±0.09.ys. 


4.1 ±0.5.,,., 


±0.7,y„ 


l.Il±0.03„., 


±0.07,y„ 


Superior conjunction (angle cut) 


2.15±0.03.„a, 


±0.09sys. 


5.0 ±0.7,,., 


±0.7,y„ 


0.91 ± 0.02,,., 


±0.07,y„ 


Periastron 


2.14 ±0.02.,., 


±0.09.ys. 


4.1 ±0.4,,., 


±0.7,y„ 


1.01 ±0.02,,., 


±0.07,y„ 


Apastron 


2.01 ±0.03,,., 


±0.09sy„ 


3.7 ±0.3,,., 


±0.7,y„ 


0.90 ± 0.02,,., 


±0.07,y„ 
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Fig. 6. — The Lomb-Scargle periodogram of the whole LS I +61°303 GeV 
data set. The orbital period is clearly visible. 



LS I +61°303 (lGregorvll2002l) . The lightcurve clearly shows 
long-term variability. We searched for changes in the orbital 
modulation of LS 1 +6 1 °303 by dividing the aperture photom- 
etry lightcurves into 6-month segments and calculating the 
power spectrum of each segment, as shown in the left panel of 
Figure|7] The amplitude of the orbital modulation is estimated 
by fitting a sine wave fixed to the orbital period to each of the 
lightcurve segments; the results, as seen in the right panel of 
Figure|7] clearly show a decreasing trend in the orbital modu- 
lation with time. 



LS I H-61°303 is one of the brightest sources in the 7-ray 
sky and towers above all other emitters in its neighborhood. 
This allows us to compute a lightcurve with an orbital bin- 
ning (26.496 days per bin) which is shown in Figure|5] Even 
by eye it is clear that the source is highly variable on orbital 
time scales and longer The longer term trends are evident by 
looking at a plot of the 3 -orbit rolling average (black line in 
Fig. |5]l. During the first eight orbits the flux decreases by a 
factor of ~2. Then, in March 2009, the flux appears to in- 
crease over the course of several orbits; we take the transition 
point of this increase to be MJD 54900. 

The flux increases significantly by 33±4%, rising from a 
baseline of (0.75 ± 0.03stat ± 0.07 ^yst) x lO-^'phcm-'s"' ob- 
tained from the first 8 months of data to (LOO ± O.Olstat ± 
0.07syst) X 10~^phcm"- s"' which is the average flux of the re- 
maining 1 .7 years of the data. Comparing the flux levels aver- 
aged over the same time span, 8 months before and 8 months 
after the flux change, we obtain a ^40% increase. After this 
flux change the flux decreases again slowly over the remain- 
ing L7 years. The complexities of the short timescale, orbit- 
to-orbit variability make it impossible to characterise the ex- 
act properties of the transition from the 'lower' to 'higher' 
flux states. The transition likely took place over several or- 
bits, however, for simplicity throughout the remainder of this 
analysis we use a transition time of MJD 54900. 

We graphically show the flux change in Figure [H by plot- 
ting the folded lightcurves before and after the tra nsition in 
March 2009. The data points are folded on the iGregorvl 
(I2OO2I) period, with zero phase at MJD 43,366.775. Before 
the transition, the modulation was clearly seen and is com- 
patible with the already published phasogram, whereas af- 
terwards, the amplitude of the modulation diminishes. We 
quantify this behavior by measuring the flux fraction below. 
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Fig. 7. — Left: The Lomb-Scargle periodograms of the LS I +61°303 30 months Hghtcurve split into five consecutive segments; the earUest is at the top. The 
red dashed line indicates the known orbital period. Right: The change in the amplitude and phase of the orbital flux modulation calculated by fitting a sine wave 
to each of the lightcurve segments. 
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Fig. 8. — Left: Folded lightcurve (lOOMeV-300 GeV) of LS I -f61°303 before the flux change (blue), when the modulation is still clearly visible, compared 
with the earher published 8-months dataset shown in gray. Middle: Folded Hghtcurve ( 100 MeV-300 GeV) after the flux change, in March 2009 (red). The 
modulation gets fainter. For comparison, the previous published lightcurve is also plotted in gray. Right: Comparison of the spectra derived before (blue) and 
after (red) the flux change in March 2009. 



Note that the data sets corresponding to the reported resuhs 
dAbdo et al.l2009al) and what we here referred to as before the 
flux change span almost exactly the same time range, with the 
consequence of our current analysis essentially reproducing 
that previously published. The time span covered by our ear- 
lier publication coincidentally finished just prior to the onset 
of the flux change. The spectra derived before and after this 
flux change are shown in Figure [8] where the increase in flux 
is also obviously visible. 

4.3. Phase resolved spectral analysis 

The statistics of the current dataset allow us to divide the or- 
bit in different phase ranges and to compute the corresponding 
spectra for different phase bins. We have divided the orbit into 
INFC and SUPC phase ranges in two different ways. First, we 
have split the LS I H-61°303 orbit in two halves based on its 
geometry, as visualized, for instance in Aragona et al. (2009). 
The SUPC phase range is defined as from phase 0.63 to phase 
0.13; INFC is defined correspondingly, as the remaining half. 
We also adopted another way to separate between INFC and 



SUPC based on the angle between the compact object, the 
star, and the observer. Therefore, the orbit is not divided in 
two halves in this way, but in one piece when the compact 
object is in front of the star (corresponding to INFC): 0.244 
- 0.507; and in another phase range with the same duration, 
centering at the exact SUPC phase: 0.981 - 0.244. As can 
be seen in Figure |9] the spectra obtained with the different 
cuts do not differ significantly and all spectral parameters are 
compatible within their corresponding errors (see Table |2]). 
We find that the flux difference between INFC and SUPC is 
of order of 20%. For these different datasets, representing 
only portions of the orbit, we modeled the source with a pure 
power law and with a power law with an exponential cutoff. 
The ATS value comparing both fits for the angular based cut 
is 108 (INFC), which means that the probability of incorrectly 
rejecting a power law with respect to an an exponential cutoff 
is 2.8 X 10-25 (lOcr). For SUPC the ATS value is 97, which 
leads to a probability of 6.9 x lO''^ (6.8f7) to wrongly reject 
the cutoff power law. Hence, the exponential cutoff is pre- 
ferred over the pure power law also for two parts of the orbit. 
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Fig. 9. — Spectra of LS I +61°303 for different phase bins during the orbit are shown. Left: Spectra for the geometrical cut of the orbit in superior (blue) and 
inferior (red) conjunction. Middle: Spectra for the angle based cut in superior and inferior conjunction. Right: Spectra derived for periastron (red) and apastron 
(blue). 
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Fig. 10. — The modulation of the spectral pai'amete rs of LS I +61°303 on 
the total 30 months lightcurve. For details see Section l431 



namely INFC or SUPC. 

We have also divided the orbit into phase ranges corre- 
sponding to periastron (half of the orbit around phase 0.275) 
and apastron (the other half, around phase 0.725), based just 
on the distance between the compact object and the star. In 
this case, neither a significant difference between the flux val- 
ues for the two phase bins nor a difference in the spectral 
shape is visible. These results can also be seen in Figure |9l 
This is probably the result of dividing the orbit into phase 
ranges which contain both bin phases corresponding to INFC 
and SUPC. 

We also studied the spectral behavior of the source in phase 
bins of 0.1. For this study we modeled LS I -i-61°303 with a 
power law with cutoff for each phase-bin individually. The 
spectral parameters obtained are shown in Figure [TO] Note 
that we fixed in the model the index at 2.07 to study the orbital 
behavior of the cutoff energy and we fixed the cutoff energy at 
3.9 GeV to study the behavior of the index, since both param- 



eters are correlated. These values are, respectively, the results 
arising from the fit over the whole dataset. A clear orbital 
modulation of the flux and spectral shape are seen. Through 
the periastron passage the spectrum gets softer and the flux 
is maximum, whereas around apastron the spectrum becomes 
harder and the flux reaches its minimum. 

4.4. Spectral fitting 

Both the orbitally-averaged spectrum and the spectra of the 
several datasets mentioned were also fit with a broken power 
law, in addition to the pure and exponentially cutoff power 
law. All the TS values for the different fits are listed in Ta- 
ble[T] It is evident that the power law with an exponential cut- 
off or the broken power law always describe the spectral shape 
better than a pure power law does. To be precise, when com- 
paring the exponentially cutoff with the pure power law, the 
ATS values span the range from 90 to 520; with the former 
being always statistically preferred. Fitting a broken power 
law gives almost the same results. The ATS values span the 
range from 77 to 494 and we find break energies in the range 
of 0.4 to 1.7 GeV. All of the fits to the various data sets with 
an exponentially cutoff power law have slightly better TS val- 
ues with fewer degrees of freedom than a broken power law 
suggesting that the former model is a better description of the 
data than the latter one. The exponentially cutoff power law 
describes nicely the curvature of the spectrum especially at 
low energies. At higher energies, above the spectral break, the 
broken power law fits all the spectral points, even the highest 
one which might be considered possibly part of a second com- 
ponent. Statistically we cannot distinguish which of these two 
fitting models describes best the data as a whole, but only that 
both of them are preferred over a pure power law. 

4.5. The multi-wavelength context 
4.5.1. X-rays 

LS I H-61°303 has also been monitored with the RXTE- 
Proportional Counter Array (PCA) and folded lightcurves 
wer e pro duced using the same ephemeris as described in Sec- 
tion l4.2l In the left panel of Figure[TTl we show a direct com- 
parison between the phasograms in X-ray and in 7-rays, with 
simultaneously taken data. We divide the whole LAT dataset 
into five periods of six months each and compare them with 
correspondingly obtained PCA X-ray data. The division in 
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Fig. 1 1 . — Comparison between the 7-ray (blue) and the X-ray (gray) data of LS I -1-61° 303 . Left: For each of the 5 separate 6-months periods, the 3-10 keV 
(gray) and the 100MeV-300GeV (blue) folded lightcurves are shown. Right: 6-months modulated flux fraction data in both energy bands, with equal color 
coding. 



periods of 6-months is justified in order to have enough statis- 
tics in 7-rays for each individual time bin, and such that orbit- 
to-orbit X-ray variability does not dominate the flux fraction 
changes. 

For RXTE-PCA, we have used the "Standard 2" mode for 
spectral analysis. Data reduction was performed using HEA- 
SOFT 6.9. We select time intervals where the source elevation 
above Earth limb is > 10° and the pointing offset is < 0.02°. 
PCA background lightcurves and spectra were generated us- 
ing the FTOOLS task pcabackest. pcarsp was used 
to generate PCA response matrices for spectra. The back- 
ground file used in analysis of PCA data is the most recent 
available from the HEASARC website for faint sources, and 
detector breakdown events have been removecQ. A power 
law shape, with absorbing hydrogen column density fixed 
at 0.75 X 10-2 cm-2 (Kalberla etaI1l2005l: iSmith et al.ll2009h . 
was used to fit the Standard 2 data, with the following function 
N{E) = Ke-^"''^'^\E/keW)-", where K is a normalization at 1 
keV, a is the photoelectric cross section and a is the photon 
index. 

It is apparent that the X-ray modulation is always visi- 

The background file is pca_bkgd_cmfaintl7_eMv20051128.mdl an d see 
the website: http://heasarc.gsfc.nasa.gov/docs/xte/recipes/pca_breakdownl itmll 
for more information on the breakdowns. The data have been barycentereH 
using the FTOOLS routine faxbary using the JPL DE405 solar system 
ephemeris. 



ble in each of these five panels, albeit with variable am- 
plitude of flux modulation. We define the flux fraction as 
(cmax-Cnvn)/(cmax + c,„in), whcrc c„,„ and c„„„ are the maxi- 
mum and minimum flux in the 3-10 keV found in the orbital 
profile analyzed (after background subtraction). The modu- 
lation becomes even stronger over time and stays stable over 
the last 3 half year bins. Instead, at GeV energies, the LAT 
data indicate that the modulated fraction fades away until the 
variability along the orbit is barely visible in the last 6 months 
of our data, which is consistent with Figure [T] The flux frac- 
tion is plotted in the right-top panel of Figure [TT] The differ- 
ence between the behavior of the X-ray and the 7-ray emis- 
sion is clearly visible. In the right-bottom panel of Figure [TT] 
we show the average 7- and X-ray fluxes fitted in each of the 
periods considered. We checked for a possible appearance of 
the super-orbital periodicity 1667±8 days, taking zero phase 
at MJD 43366.275 and the shape of the outburst peak flux 
modulation estimated by Gregory (2002) in radio by consid- 
ering the flux evolution and the direct count rate, using the full 
dataset of the PCA observations. We find no evidence of the 
super-orbital period in X-ray or 7-rays, which is not surpris- 
ing given the short integration time in comparison with the 
super-orbital period duration. 

4.5.2. Radio and optical 
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Fig. 12. — Comparison of non-contemporaneous LAT-obtained 7-ray with radio data of LS I +61°303 folded with the radio super orbital period. Top: With 
respect to the radio data compiled by Gregory (1999), mostly at 8.3 GHz. Bottom: With respect to the EW of the Ha line obtained by Zamanov et al. (1999). 



We have also compared the LAT-obtained data with radio 
and Ha observations. For the latter, we first take into account 
the results contained in the long-term coverage presented by 
IZamanov et al.l ( ll999l) . These authors have already shown that 
the equivalent width (EW) and the peak separation of the Ha 
emission line appear to vary with the super-orbital radio pe- 
riod of ~ 1600 days, this being likely the result of cyclical 
variations in the mass-loss rate of the Be co mpanion and/or 
of density variability in the circumstellar disk. IZamanov et alj 
{1999) proposed that the variability in the EW(Ha) can be 
explained by a cyclical change in the mass-loss rate of about 
25% over its average value. This mechanism would imply 
changes in the density of the circumstellar disc in the same 
range. The fact that the 7-ray lightcurve would have a simi- 
lar behavior to that found for the EW(Ha) would imply that 
the former is related to the circumstellar disc surrounding the 
Be star This would naturally be the case if the 7-rays are 
produced in an inter-wind shock formed by the colli sion of 
outflow s fr om the compact object and the s tar its elf. iDubusI 
(120061) and iSierpowska-Bartosik & TorresI (120091) discussed 
how the many uncertainties present in our knowledge of the 
LS 1 H-61°303 system influence the models for producing 7- 
ray emission. One of the important parameters is the density 
and size (and the possible truncation) of the circumstellar disc. 
The latter influence the position of the inter-wind shock, and 
the time intervals along the orbit in which the compact ob- 
ject outflow may be balanced by the equatorial wind feeding 
it, and/or when the compact object is directly within the cir- 
cumstellar disc itself. For recent measurements on very long- 
term optical variability of Be High Mass X-ray Binaries i n the 
Small Magellanic Cloud see lRaioelimanana et al.l (1201 Ih . 

We plot in Figure [12] a comparison of the LAT data, folded 
on the super-orbital radio period, with the radio data com- 
piled by Gregory (1999) - most of it obtained at 8.3 GHz, 



using the Green Bank Interferometer - and the EW of the 
Ha line obtained by IZamanov et al.l ( 119991) . Care should 
be exer cised when comparing our LAT data with the radio 
data in IZamanov et al.l d 19991) s ince the latter authors used 
the ephemeris given by iGregorvl (11999), for which the com- 
bination of the phase and peak flux density yielded a best- 
modulation period o f 1584 d a ys. T his super-orbital period 
was later revised by IGregorvl (120021) to the current value of 
1667 days. This latter value is the super-orbital period that 
we have used to fold the X-ray and LAT data, and we use it 
also here for the radio and Ha emission. The comparison is 
clearly non-contemporaneous, and this may induce problems 
in the interpretation of a source like LS 1 H-61°303, presenting 
different variability phenomenology. We do not see any clear 
correlation or anti-correlation with the radio and Ha fluxes 
which may be hidden by the scarcity of GeV data. 

We have also compared the Fermi data with simultaneous 
radio data from the Owens Valley Radio Observatory (OVRO) 
and the Arcminute Microkelvin Imager (AMI) array (Cam- 
bridge, UK). Since the launch of Fermi in June 2008, the 
OVRO 40 m single-dish telescope, that is located in Califor- 
nia (USA), has conducted a regu lar monitoring prog ram of 
Galactic binaries (e.g. Cyg X-3, lAbdo et alj (12009c')). The 
OVRO flux densities are measured in a single 3 GHz wide 
band centered on 15 GHz. A complete description of the 
O VRO 40 m telescope a nd calibration strategy can be found 
in iRichards et al.l ( 1201 ll) . Furthermore, we also use comple- 
mentary observations provided by AMI, consisting of a set of 
eight 13 m antennas with a maximum baseline of ~ 120 m. 
AMI observations are conducted with a 6 GHz bandwidth re- 
ceivers also centered at 15 GHz. See iZwart et al.l (120081) for 
more details on the AMI interferometer, that is mostly used 
for study of the cosmic microwave background. By folding 
these radio data no super-orbital modulation could be seen. 
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Fig. 13. — Lightcurve of LS I +61° 303. Each point represent one orbit of 26.496 days. The black dashed line marks the moment of the flux change in March 
2009, together with the simultaneously taken radio data by the AMI and OVRO instruments (both operating at 15 GHz) are shown in gray. 



which could be due to the poor coverage of a whole super- 
orbit. A direct c omp arison of the GeV and the radio data is 
shown in Figure [T3l where the flux over time is plotted. No 
correlation of the data points can be seen either In summary, 
we did not find any correlation between the GeV and the radio 
band, either in archived nor in simultaneous data. 

5. CONCLUDING REMARKS 

After analyzing a dataset comprising 2.5-years of Fermi- 
LAT observations of the two binaries LS 5039 and 
LS 1 +61°303, we note several cha nges with respect to the 
initial reports of lAbdo et alj (I2009allbl) . These were produced 
either because the accumulation of a longer observation time 
allowed us to make distinctions that were earlier impossible 
(valid for both sources), or because the behavior of the source 
changed (valid for LS 1 +61°303). On one hand, the statis- 
tics are now sufficient to divide the dataset of both sources in 
INFC and SUPC and to show that a power law with a cutoff 
describes the spectra obtained in both conjunctions better than 
a pure power law. The cutoff is similar to that found in the 
many other GeV pulsars discovered by fermZ-LAT. However, 
we have found that both LS 5039 and LS 1 H-61°303 show an 
excess of the high energy GeV emission beyond what is ex- 
pected from an exponentially cutoff power law. While the 
high energy data are significantly in excess of the exponen- 
tially cutoff power law there are insufficient statistics at these 
high energies to model the excess with an additional spectral 
component. 

The process(es) generating such a component in the case 
of LS5039 and LS 1 H-61°303 is unclear and may even be 
different in the two sources. However, such a second com- 
ponent would present a possible connection between the GeV 
and TeV spectra in both sources. Collecting more data and 
therefore more statistics will allow to prove or discard it in 



the future. The lack of datapoints at high energies, also af- 
fects, particularly in the case of LS 1 H-61°303, the distinction 
between an exponentially cut and a broken power law. Cur- 
rently, both are certainly preferred over a pure power law, but 
differences in the significances provided between the former 
are minor 

We have noted that whereas LS 5039 shows stable emission 
over time and also a stable orbital modulation, LS 1 H-61°303 
shows a change in flux in March 2009. Afterwards the orbital 
modulation decreases (see the bottom-most left panels of Fig- 
urefTTT) and the orbital period could not be detected in the GeV 
data. LS 1 H-61°303 has also presented a complex, concurrent 
behavior at higher energies. At TeV, for approximately the last 
two years, it seems to have been in a low state in comparison 
with the flux level t hat led to its discovery ( lAleksic et al.l201 it 
lAcciari et alj|201 11) . Additionally, it was detected once -after 
4.2 hrs of observations- near the periastr on, where the system 
was never seen at TeV energies before (lAcciari et al.ll201 lb . 
Both of these aspects of the LS 1 H-61°303 phenomenology - 
as well as the phase location of the TeV maxima in general- 
would require modifications of simple inverse Compton mod- 
els in the scenarios usually put forward for the source. The 
idea of a magnetar co mpact object in LS 1 H-61°303 is sug- 
gested bv lTorres et al.i ( 120121) to explain the changing TeV be- 
havior of this object and could potentially describe the dimin- 
ishing orbital modula tion now observed a t GeV fluxes. How- 
ever, as explained by [Torres et"al] ( 120121) . detailed numerical 
simulations are required to verify if this scenario can accu- 
rately reproduce the observed GeV and TeV emission. 

The lower-energy, multi-wavelength picture of 
LS 1 H-61°303 is still unclear with the orbital modula- 
tion in GeV (X-ray) fading (increasing) with time. Because 
the GeV data do not cover a whole super-orbital period, 
a conclusion about the possible relationship of the GeV 
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emission with the super-orbital behavior of LS I +61°303 
cannot be drawn yet. Continued monitoring of the source 
with the Fermi-LAiY will allow the study such long-term 
cycles, if any. Zamanov et al. (1999) concluded that if the 
^4yr modulation in radio and Ha is due to changes in the 
circumstellar disk density, it could be detected in X- or 7-rays 
as well, since the high energy emission of a putative neutron 
star depends on the density of the surrounding matter with 
which it interacts. Indeed, thi s modulation was recently 
revealed in X-rays (ILi et alJI2012l) . 

The nature of these systems is thus still unclear and multi- 
wavelength observations, including ferm/-LAT monitoring, 
as well as further theoretical studies are expected to bring new 
insights. 
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